Metaprogramming: "Code that creates code"

Julia has strong metaprogramming capabilities, i.e. we can use Julia to manipulate objects in Julia that represent Julia code. In this way, we can produce code in an automatic way.

Since this is rather abstract, let's consider a very concrete example: Wilkinson-type polynomials. The Wilkinson polynomial is

$$p_{20}(x) := (x-1) \cdot (x-2) \cdot \cdots \cdot (x-20) = \prod_{i=1}^{20} (x-i)$$

[Polynomials like this are interesting since their eigenvalues are very sensitive to perturbations of the coefficients of the polynomial.]

Suppose we wish to define this polynomial in Julia. The simple way would be to write it explicitly:



In [ ]:

    
p_5(x) = (x-1)*(x-2)*(x-3)*(x-4)*(x-5)

$p_{10}$ is already a pain to type by hand, $p_{20}$ more so, and $p_{100}$ is basically impossible. But this is just a case of repetition, and computers are designed for that. Of course, one possible definition uses a for loop:



In [ ]:

    
function wilkinson(n, x)
    result = x - 1
    
    for i in 2:n
        result *= (x - i)
    end
    
    result
end

We can use an anonymous function to have the actual function object $p_n$:



In [ ]:

    
wilkinson(n) = x -> wilkinson(n, x)

However, as we have seen, currently "anonymous functions are slow". It seems like it should be possible to use the original definition of the function by "unrolling the loop" and writing a loop that writes the code to generate the function.

Expressions

In other languages, e.g. Python, code is manipulated as strings. Julia takes a different approach. Consider the string



In [ ]:

    
s = "(x-1) * (x-2)"

To convert this into an object that represents Julia code, we parse it:



In [ ]:

    
ex = parse(s)



In [ ]:

    
typeof(ans)

ex is a Julia expression object. It can be thought of as a representation of the "abstract syntax tree" (AST) representing the internal structure of the expression. We can see this in two ways, using dump:



In [ ]:

    
dump(ex)

which shows everything in detail, or



In [ ]:

    
Meta.show_sexpr(ex)

which gives a compact version.

We see that Exprs have a hierarchical format that represents the code. Since they are simply Julia objects, e.g. Arrays, we can use Julia to manipulate them.

The object :x is of type Symbol, and represents the unevaluated object :x. This is basically a special type of string.

Exercise: Write a function that takes an expression object and replaces all the :xs by :zs.

Code interpolation

In our case, we do not need to mess with the internal structure of code, but rather build up code from preexisting bits. We can do this as follows:



In [ ]:

    
ex = :(x-1)



In [ ]:

    
ex = :(($ex) * (x-2))

Here, in the same way as in string interpolation, we have inserted the current value of ex into the expression!



In [ ]:

    
ex = :(($ex) * (x-3))

Now we can make our loop:



In [ ]:

    
n = 10
ex = :(x-1)
for i in 2:n
    ex = :( ($ex) * (x-i) )
end



In [ ]:

    
ex

This did not work, since we did not want "the code 'i'", but rather the value of i. So:



In [ ]:

    
n = 10
ex = :(x-1)
for i in 2:n
    ex = :( ($ex) * (x-$i) )
end



In [ ]:

    
ex

Now we need to produce the name of the function:



In [ ]:

    
name = symbol(string("W_", n))

[In Julia v0.4, this can be written more simply as symbol("W_", n).]

So the code that we would like is



In [ ]:

    
code = :( $name(x) = $ex )

Now we wish to evaluate this:



In [ ]:

    
eval(code)

This creates a function with the name W_10 that does (almost) exactly what we would write by hand.

Let's compare the two options by evaluating the function on a grid of values:



In [ ]:

    
f1(range) = [W_10(x) for x in range]
f2(range) = [wilkinson(10, x) for x in range]



In [ ]:

    
range = -10:0.000001:10
@time f1(range);
@time f2(range);

We see that the generated code with the unrolled loop is twice as fast as the naive for loop.

Code generation like this is used frequently in Julia code when repetitive code is required.

Macros

The things that feel like functions whose names start with @ that we have been using, such as @time, @which, etc., are not functions in the standard sense, but rather are macros. These are "super-functions" whose arguments are code expressions, and which transform one code expression into another. This new code expression is then evaluated as if the new code had been typed directly.

To understand what is going on, let's define a simple macro:



In [ ]:

    
macro simple(expr)
    @show expr
    0  # for the moment
end



In [ ]:

    
@simple x+y

We see that the Julia code that follows the macro call is passed to the macro already having been parsed into an Expr object.

Suppose we redefine the macro as



In [ ]:

    
macro simple(expr)
    @show expr
    expr  # returns expr
end

Then we get



In [ ]:

    
@simple x+y

What is happening here is that the macro returns the expression :(x+y), and this is then evaluated using eval. The result is that Julia tries to calculate the value of the expression x+y, but the variable x is not defined. Let's define y and z, but not x:



In [ ]:

    
y = 3; z = 4



In [ ]:

    
x



In [ ]:

    
@simple x+y

Now let's define a new macro @replace that uses our previous replace function:



In [ ]:

    
macro replace!(expr)
    replace(expr)
    @show expr
    expr
end



In [ ]:

    
@replace! x+y

`macroexpand`

We can discover what a macro does using the macroexpand function which takes a code expression:



In [11]:

    
macroexpand(:(@time sin(10)))









    Out[11]:





quote  # util.jl, line 53:
    local #464#b0 = Base.gc_bytes() # line 54:
    local #465#t0 = Base.time_ns() # line 55:
    local #466#g0 = Base.gc_time_ns() # line 56:
    local #467#val = sin(10) # line 57:
    local #468#g1 = Base.gc_time_ns() # line 58:
    local #469#t1 = Base.time_ns() # line 59:
    local #470#b1 = Base.gc_bytes() # line 60:
    Base.time_print(Base.-(#469#t1,#465#t0),Base.-(#470#b1,#464#b0),Base.-(#468#g1,#466#g0)) # line 61:
    #467#val
end

As an example of the use of a macro, let's look at the @interval macro from the ValidatedNumerics.jl package:



In [ ]:

    
Pkg.add("ValidatedNumerics")
Pkg.checkout("ValidatedNumerics", "master")  # checkout the master branch ["master" is not necessary]



In [12]:

    
using ValidatedNumerics



In [13]:

    
@interval(0.1)









    Out[13]:





[0.09999999999999999, 0.1]

Floating-point calculations are not precise. Interval arithmetic is a way to provide a guaranteed enclosure of a result, i.e. an interval that contains the true result.



In [14]:

    
macroexpand(:(@interval(0.1)))









    Out[14]:





:(ValidatedNumerics.make_interval(interval_parameters.precision_type,0.1))



In [15]:

    
@interval sin(0.1) + cos(0.2)









    Out[15]:





[1.0798999944880696, 1.07989999448807]



In [16]:

    
macroexpand(:(@interval sin(0.1) + cos(0.2)))









    Out[16]:





:(ValidatedNumerics.+(sin(ValidatedNumerics.make_interval(interval_parameters.precision_type,0.1)),cos(ValidatedNumerics.make_interval(interval_parameters.precision_type,0.2))))

This has the structure sin(interval(0.1)) + cos(interval(0.2)) (although the details are a bit more complicated).

So basically the @interval macro does something very similar to our replace! macro: it walks the tree of the expression; finds parts of it that are of a certain type (numeric literals); and wraps those in the make_interval function